This class provides an interface both to the ASE environment and the Fortran core of Pysic through the pysic_fortran module. In ASE terms, Pysic is a calculator, i.e., an object that calculates forces and energies of given atomistic structures provided as an ASE Atoms object.
Pysic is not a fixed potential calculator, where the interactions are determined solely based on the elements present in the system. Instead, Pysic allows and requires the user to specify the interactions governing the system. This is done by providing the calculator with one or several Potential objects.
Calculators in ASE are a class representing the physics of the atomistic system. Or put in another way, a calculator is the mathematical machine containing the algorithms and routines for calculating, for instance, the total potential energy of the atomistic system represented as an ASE atoms object (hence the name). As ASE itself contains only some very simple calculators, in practice the calculators are often just Python interfaces to other atomistic simulators. Pysic too is such a calculator.
Pysic implements classical potentials, where the interactions between atomic pairs, triplets etc. are represented via simple mathematical functions or a Tabulated potential. These can be combined freely using the Potential class and its modifiers and extensions such as the Coordinator, Bond order parameter, Product potential, Compound Potential, and Coulomb summation classes.
Although Pysic can be used as a calculator on its own, it is also possible to use it as a wrapper for other ASE calculators. By doing so, one may combine different simulators and even different descriptions in one calculator, and all the wrapped calculators will be executed when the Pysic calculator needs to calculate energies or forces. This is done simply by adding the external ASE calculators to Pysic using the add_calculator() method.
As a simple example, one can use Pysic to add simple components like springs in an ab initio calculation. This example uses the Gpaw code.:
import pysic
from ase import Atoms
from gpaw import GPAW
# create the system
system = Atoms( ... )
# create the pysic calculator
calc = pysic.Pysic()
# create a classical spring potential between atoms '0' and '1'
spr = pysic.Potential('spring', indices=[0,1], parameters=[10.0, 3.0], cutoff=10.0, cutoff_margin=1.0)
calc.add_potential(spr)
# create an ab initio calculator
dft = GPAW(h=0.18, nbands=1, xc='PBE', txt='wrapped.out')
# wrap the calculators
calc.add_calculator(dft)
# calculate the energy
system.set_calculator(calc)
print system.get_potential_energy()
It is also possible to wrap several instances of Pysic together, but usually there should be no need to do this. It is more efficient to just add all the potentials in one Pysic instance. However, as an example, here is how it works:
import pysic
import ase
# create some test system
system = ase.Atoms('COH',[[0.0,0.0,0.0], [2.0,0.0,0.0], [2.0,2.0,2.0]])
# create two calculators
calc1 = pysic.Pysic()
calc2 = pysic.Pysic()
# potential 1
pot1 = pysic.Potential('LJ',symbols=['C','O'],parameters=[1,1],cutoff=10.0)
# potential 2
pot2 = pysic.Potential('LJ',symbols=['O','H'],parameters=[1,1],cutoff=10.0)
# add potentials to respective calculators
calc1.add_potential(pot1)
calc2.add_potential(pot2)
system.set_calculator(calc1)
# write energies and forces for the calculators both separately and combined
print "Only potential 1"
print system.get_potential_energy()
print system.get_forces()
print "\n"
print "Only potential 2"
print calc2.get_potential_energy(system)
print calc2.get_forces(system)
print "\n"
print "Both potentials by calculator wrapping"
calc1.add_calculator(calc2)
print system.get_potential_energy()
print system.get_forces()
print "\n"
calc1.remove_calculator(calc2)
calc1.add_potential(pot2)
print "Both potentials in one calculator"
print system.get_potential_energy()
print system.get_forces()
print "\n"
This will output:
Only potential 1
-0.015380859375
[[ 0.04541016 0. 0. ]
[-0.04541016 0. 0. ]
[ 0. 0. 0. ]]
Only potential 2
-0.00194931030273
[[ 0. 0. 0. ]
[ 0. 0.00291824 0.00291824]
[ 0. -0.00291824 -0.00291824]]
Both potentials by calculator wrapping
-0.0173301696777
[[ 0.04541016 0. 0. ]
[-0.04541016 0.00291824 0.00291824]
[ 0. -0.00291824 -0.00291824]]
Both potentials in one calculator
-0.0173301696777
[[ 0.04541016 0. 0. ]
[-0.04541016 0.00291824 0.00291824]
[ 0. -0.00291824 -0.00291824]]
Notice how in the third printout the calculators calc1 and calc2 are wrapped together, and the result is their sum. In the fourth case, calc2 is discarded and instead pot2 is added directly to calc1. The result is the same, but the latter option is better performance wise, since wrapping the calculators requires both to be separately managed during calculation.
Below is a list of methods in Pysic, grouped according to the type of functionality.
A calculator class providing the necessary methods for interfacing with ASE.
Pysic is a calculator for evaluating energies and forces for given atomic structures according to the given Potential set. Neither the geometry nor the potentials have to be specified upon creating the calculator, as they can be specified or changed later. They are necessary for actual calculation, of course.
Simulation geometries must be defined as ASE Atoms. This object contains both the atomistic coordinates and supercell parameters.
Potentials must be defined as a list of Potential objects. The total potential of the system is then the sum of the individual potentials.
Parameters:
Add a calculator to the list of external calculators.
Parameters:
Add a potential to the list of potentials.
Also a list of potentials can be given as an argument. In that case, the potentials are added one by one.
If a CompoundPotential is given, the addition is done through its build() method.
Parameters:
Calculates electronegativities.
Calls the Fortran core to calculate forces for the currently assigned structure.
Calculates the potential energy.
Calls the Fortran core to calculate the potential energy for the currently assigned structure.
If a link exists to a ChargeRelaxation, it is first made to relax the atomic charges before the forces are calculated.
Calculates forces (and the potential part of the stress tensor).
Calls the Fortran core to calculate forces for the currently assigned structure.
If a link exists to a ChargeRelaxation, it is first made to relax the atomic charges before the forces are calculated.
Calculates the potential part of the stress tensor (and forces).
Calls the Fortran core to calculate the stress tensor for the currently assigned structure.
Check if a calculation is required.
When forces or energy are calculated, the calculator saves the result in case it is needed several times. This method tells if a wanted quantity is not yet calculated for the current structure and needs to be calculated explicitly. If a list of several quantities is given, the method returns true if any one of them needs to be calculated.
Parameters:
An object storing the data passed to the core.
Whenever a Pysic calculator alters the Fortran core, it should also modify the core object so that it is always a valid representation of the actual core. Then, whenever Pysic needs to check if the representation in the core is up to date, it only needs to compare against core instead of accessing the Fortran core itself.
Returns true if the core is always fully initialized, false otherwise.
Initializes the neighbor lists.
In order to do calculations at reasonable speed, the calculator needs a list of neighbors for each atom. For this purpose, either the list is built in the Fortran core or the ASE NeighborList are used. The ASE lists are very slow, but they will work even for small systems where the cutoff is longer than cell size. The custom list uses pysic.calculator.FastNeighborList, but the build succeeds only for cutoffs shorter than cell size. The type of list is determined automatically.
This method initializes these lists according to the given cutoffs.
Parameters:
Set the core initialization mode.
Parameters:
Returns the list of other calculators to be used with Pysic.
Returns the ChargeRelaxation object connected to the calculator.
Update for ASE 3.7
Returns the electronegativity differences of atoms from the average of the entire system.
Returns the forces.
If the atoms parameter is given, it will be used for updating the structure assigned to the calculator prior to calculating the forces. Otherwise the structure already associated with the calculator is used.
The calculator checks if the forces have been calculated already via calculation_required(). If the structure has changed, the forces are calculated using calculate_forces()
Parameters:
Get a list of maximum cutoffs for all atoms.
For each atom, the interaction with the longest cutoff is found and the associated maximum cutoffs are returned as a list. In case the a list of scaled values are required, the scaler can be adjusted. E.g., scaler = 0.5 will return the cutoffs halved.
Parameters:
Returns the neighbor list object.
Returns the FastNeighborList or ASE NeighborList object assigned to the calculator.
The neighbor lists are generated according to the given ASE Atoms object and the Potential objects of the calculator. Note that the lists are created when the core is set or if the method create_neighbor_lists() is called.
Numerically calculates the gradient of a bond order factor with respect to moving a single particle.
This is for debugging the bond orders.
Numerically calculates the derivative of energy with respect to charging a single particle.
This is for debugging the electronegativities.
Numerically calculates the negative gradient of energy with respect to moving a single particle.
This is for debugging the forces.
Returns the potential energy.
If the atoms parameter is given, it will be used for updating the structure assigned to the calculator prior to calculating the energy. Otherwise the structure already associated with the calculator is used.
The calculator checks if the energy has been calculated already via calculation_required(). If the structure has changed, the energy is calculated using calculate_energy()
Parameters:
Returns the stress tensor in the format \([\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{yz},\sigma_{xz},\sigma_{xy}]\)
If the atoms parameter is given, it will be used for updating the structure assigned to the calculator prior to calculating the stress. Otherwise the structure already associated with the calculator is used.
The calculator checks if the stress has been calculated already via calculation_required(). If the structure has changed, the stress is calculated using calculate_stress()
Stress (potential part) and force are evaluated in tandem. Therefore, invoking the evaluation of one automatically leads to the evaluation of the other. Thus, if you have just evaluated the forces, the stress will already be known.
This is because the stress tensor is formally defined as
where \(m\), \(v\), \(r\), and \(f\) are mass, velocity, position and force of atom \(i\), and \(A\), \(B\) denote the cartesian coordinates \(x,y,z\). (The minus sign is there just to be consistent with the NPT routines in ASE.) However, if periodic boundaries are used, the absolute coordinates cannot be used (there would be discontinuities at the boundaries of the simulation cell). Instead, the potential energy terms \((r_i)_A (f_i)_B\) must be evaluated locally for pair, triplet, and many body forces using the relative coordinates of the particles involved in the local interactions. These coordinates are only available during the actual force evaluation when the local interactions are looped over. Thus, calculating the stress requires doing the full force evaluation cycle. On the other hand, calculating the stress is not a great effort compared to the force evaluation, so it is convenient to evaluate the stress always when the forces are evaluated.
Parameters:
Fully initializes the Fortran core, creating the atoms, supercell, potentials, and neighbor lists.
Check if the cutoffs have been expanded.
If the cutoffs have been made longer than before, the neighbor lists have to be recalculated. This method checks the individual cutoffs of all atoms to check if the cutoffs have changed.
Parameters:
Remove a calculator from the list of calculators.
Parameters:
Remove a potential from the list of potentials.
If a CompoundPotential is given, the removal is done through its remove() method.
Parameters:
Assigns the calculator with the given structure.
This method is always called when any method is given the atomic structure as an argument. If the argument is missing or None, nothing is done. Otherwise a copy of the given structure is saved (according to the instructions in ASE API.)
If a structure is already in memory and it is different to the given one (as compared with __ne__), it is noted that all quantities are unknown for the new system. If the structure is the same as the one already known, nothing is done. This is because if one wants to access the energy of forces of the same system repeatedly, it is unnecessary to always calculate them from scratch. Therefore the calculator saves the computed values along with a flag stating that the values have been computed.
Parameters:
Add a charge relaxation algorithm to the calculator.
If a charge relaxation scheme has been added to the Pysic calculator, it will be automatically asked to do the charge relaxation before the calculation of energies or forces via charge_relaxation().
It is also possible to pass the Pysic calculator to the ChargeRelaxation algorithm without creating the opposite link using set_calculator(). In that case, the calculator does not automatically relax the charges, but the user can manually trigger the relaxation with charge_relaxation().
If you wish to remove automatic charge relaxation, just call this method again with None as argument.
Parameters:
Sets up the Fortran core for calculation.
If the core is not initialized, if the number of atoms has changed, or if full initialization is forced, the core is initialized from scratch. Otherwise, only the atomic coordinates and momenta are updated. Potentials, neighbor lists etc. are also updated if they have been edited.
Set the Coulomb summation algorithm for the calculator.
If a Coulomb summation algorithm is set, the Coulomb interactions between all charged atoms are evaluated automatically during energy and force evaluation. If not, the charges do not directly interact.
Parameters:
Copy and save the list of individual cutoff radii.
Parameters:
cutoffs: list of doubles new cutoffs
Assign a list of potentials to the calculator.
Also a single potential object can be given, instead of a list. Note that this method does not conserve any potentials that were already known by the calculator. To add potentials to the list of known potentials, use add_potential().
Parameters:
Updates the positions and momenta of atoms in the Fortran core.
The core must be initialized and the number of atoms must match. Upon the update, it is automatically checked if the neighbor lists should be updated as well.
Updates the neighbor lists in the Fortran core.
If uninitialized, the lists are created first via create_neighbor_lists().